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Summary. We investigate the effect of fragmentation on the homogeneous free 
cooling of inelastic hard spheres, using Boltzmann kinetic theory and Direct Monte 
Carlo simulations. We analyze in detail a model where dissipative collisions may 
subsequently lead to a break-up of the grains. With a given probability, two off- 
springs are then created from one of the two colliding partners, with conservation of 
mass, momentum and kinetic energy. We observe a scaling regime characterized by 
a single collisional average, that quantifies the deviations from Gaussian behaviour 
for the joint size and velocity distribution function. We also discuss the possibility 
of a catastrophe whereby the number of particles diverges in a finite time. This 
phenomenon appears correlated to a "shattering" transition marked by a delta 
singularity at vanishingly small grains for the rescaled size distribution. 



1 Introduction 

The phenomenon of fragmentation is central to many natural and technological 
processes, from geology to the pharmaceutical industry (where e.g. the properties 
of an active principle depend on the size distribution of its constituents). In the 
realm of granular gases where the fluidization is achieved by a violent excitation 
or shaking, the fragmentation mechanism is also relevant although experimental 
studies of these systems usually avoid such a complication. 

A theoretical description of fragmentation is nevertheless a difficult task and 
traditional approaches rely on rate (Smoluchowski-like) equations 0, that describe 
the time evolution of the size distribution after having integrated out the other de- 
grees of freedom, such as the velocities of the particles. However, when the motion 
between the collision events is free, the dynamics of a reacting system in general 
(where, because of collisions, particles may break up, aggregate, annihilate . . . ) is 
determined by details of the collision (collisional correlations which often in- 
validate the rate equation approach, at least in its usual version. In particular, for 
dynamics that do not conserve the number of particles, the time evolution of the 
density crucially depends on the ratio of the kinetic energy dissipated in a "typi- 
cal" collision, to the mean kinetic energy of the particles |2I3| . It therefore appears 
that a more microscopic description than the coarse grained rate equation route is 
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necessary, at the expense of a non negligible complication. When particles disap- 
pear upon colliding, it is however possible to prove that the Boltzmann equation 
becomes exact at late times (in space dimensions D > 1 |[j]). Although it does 
not seem possible to show that the validity of Boltzmann's picture (which is less 
restrictive than the Smoluchowski approach) pertains to a modification of the re- 
action rule, such an approach captures important collisional correlations missed by 
the naive rate equation framework. In the following analysis, we will therefore use 
Boltzmann's kinetic theory (relying on the molecular chaos assumption that the 2 
body distribution function of colliding pairs may be factorized in terms of the single 
particle distribution 4 ) , to investigate the dynamics of a fragmenting granular gas. 

2 Model 

We consider a gas of grains that evolve freely in time, from an initial equilibrium 
distribution, i.e. we assume that up to the initial time, the system is elastic, and 
hence, that its velocity distribution is Maxwellian. 

In the fluidized regime we shall be interested in, grains interact as hard bodies; 
collisions are always binary, and instantaneous. In a D-dimensional system, at each 
collision a fraction (1 — a 2 )/ (2D) of the relative kinetic energy is lost, where a, 
the inelasticity (or restitution) coefficient, measures the departure from elastic be- 
haviour. After each collision, the stored kinetic energy in internal grain excitations 
may lead to a break-up of the grains. Such a fragmentation event is stochastic in 
nature. Previous studies suggest than in a number of materials, the probability that 
a grain breaks is related to the number of surface defects (or micro-cracks). The 
more kinetic energy is dissipated, the larger the chance that one of such defects 
enlarges, leading to the break-up of the grain. Hence, the fracture probability is re- 
lated both to the energy lost and to the defect distribution (and the related surface 
stress distribution). There exist different models in the literature that account for 
such effects. A reasonable fragmentation probability model assumes that the stored 
energy and the defect distribution are uncorrelated, and that defects are uniformly 
distributed over a particle's surface; the breaking probability is thus proportional to 
the particle's surface 6 . Once a grain breaks, we should also specify the number of 
out-coming particles, and how their sizes and velocities are related to the physical 
properties of the fragmenting grain. Here the situation is more controversial, and 
it is highly dependent on whether the breaking process is induced by an energy 
input from a source which drives the system or whether it is due to the internal 
dynamics of the system; despite this, the Roslin-Rammler law (and a few variants) 
is considered to be a reasonable description for the off-spring distribution |H6| . 

Rather than attempting a detailed description of the fragmentation mechanism, 
in this paper we will assume the simplest break-up process. The idea is to minimize 
the amount of input at the collision level and to analyze which is the behavior of 
the system under the simplest possible mechanical rules. We hope to be able to 
carry out a detailed analysis that will be helpful before investigating more refined 
approaches. 

In our model, after the collision has taken place, with a probability p we break 
one of the two colliding grains. The grain fragments in 2 off-springs which keep the 
(spherical) shape of the parent grain. The fragmentation process conserves mass, 
momentum and kinetic energy. The mass is distributed randomly between the two 
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out-coming grains with a uniform probability density. In this case, the conservation 
of momentum and kinetic energy imply that both off-springs have the same velocity, 
equal to that of the fragmenting grain. If we would assume that the two colliding 
grains may break up due to the collision, there would exist more freedom in the 
velocity redistribution of the off-springs, even assuming coplanarity of the out- 
coming grains' trajectories. We will consider that the mass density of the grains is 
conserved, so that the relation between mass m and radius a is taken as m = a D . 

We will further concentrate on the homogeneous freely evolving state of the fluid 
under consideration. Such a situation has received much attention in the absence of 
fragmentation, and the theoretical predictions have been contrasted against molec- 
ular dynamics simulations. The homogeneous state is known to be unstable against 
low wave length fluctuations of the vorticity, and at late times, hydrodynamic flows 
set up[!j]. As a first step, we will concentrate in the regime prior to the appearance 
of the hydrodynamic instability. In fact, having a hydrodynamics origin, we might 
also expect an analogous instability in the present model. 

In this homogeneous state, the Boltzmann equation, which determines, at a 
mean-field level, the time evolution of the system through the one-particle distri- 
bution function f(v,a,t), in a D-dimensional space has the form 



In the previous expression, v stands for the velocity, a for the particles' radius, while 
e is a unit vector joining the centers of the colliding grains. The Heaviside function, 
9, ensures the appropriate kinematic constraint, that only approaching grains will 
collide, while the term \v\2 ■ e\ is the collision cross-section; the subscript 12 

refers to the difference between the values of the colliding pair, e.g. v±2 = Vi — v%- 
Moreover, g(u) is the probability density that one of the off-springs, in case of 
fragmentation, has a mass u times the parent's mass. Due to mass conservation, 
g(u) = g{l — u). In the model analyzed below, we consider the simplest case where 
this distribution is uniform [g(u) — 1 and < u < 1]. 

For clarity's sake, we have expressed the general evolution equation in terms of 
the elementary processes affecting the two colliding grains. Obviously, a simplified 
expression can be derived using the symmetric role that the two grains play during 
the collision. 

The first two terms in the r.h.s. of Q account for the disappearance of parti- 
cles of a given velocity due to the collision process itself. For the same reason, new 
particles of a given species and velocity are produced at the collision. These terms 
involve the precollision velocities, v* , because we require that the out-coming ve- 
locity of the collision corresponds to the velocity v. It is through these, precollisonal 
velocities, that the inelastic character of the collision enters. Indeed, the precollision 
velocities, in terms of the post-collisional ones (which are the relevant velocities in 
the previous kinetic equation in order to describe the kinetic processes), are given 



df(v ; a,t) = 1 



+ 



+ 



+ 




(1) 
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by 

v i,2 = «i,2 T ; (1 + a -1 ) [e ■ ui 2 ] £ (2) 

mi + m 2 v ' 

The two terms in Q containing the prefactor 1 — p correspond to the collisions in 
which, with probability 1 — p, particles do not break. The last four terms account 
for the fragmentation of grain 1 with mass mi(= erf) (resp. grain 2 with mass m 2 ), 
leading to the creation of two off springs with masses um\ and (1 — u)vn\ [resp. 
um,2 and (1 — u)m2). Mass conservation implies a correlation between the sizes of 
the newborn particles. 



3 Kinetics 

A detailed analysis of the dynamics of the system, already at the Boltzmann level, 
according to Q, is quite involved. However, the key information is already present 
in the time evolution of the total density and kinetic energy densities, defined 
respectively as 

n(t) = j f(v,a,t)dvda (3) 
E c (t) = J v 2 a D f(v,a,t)dvda. (4) 

The total mass density 

M = J a D f(v,a,t)dvda (5) 
is a conserved quantity and we define the temperature kinetically as 

2E C (mv 2 ) 

In the previous equation, the angular brackets of a quantity Q denote the mean 
value J Qf(v,a,t) dvda/n. 

From the Boltzmann equation, the time evolution for the density of particles 
reads 

^=M*Mt) (7) 

where p is the probability that a particle breaks at a collision, and ui(t) the collision 
frequency 

u)(t) = — [ dedld2 a° 2 ~ 1 9(a ■ v 12 )\i- v 12 \ f(v 1 ,ai,t)f(v 2 ,a2,t), (8) 
n J 

where the short-hand notation dld2 stands for dv\da\dv 2 da2- The kinetic energy 
density is unaffected by the breaking events (unlike the density itself), and only 
decreases due to the inelastic collisions: 

= -2b cl0 Lu{t)E c {t). (9) 

at 

In this equation, 70 = (1 — a 2 )/ (2D) is related to the energy lost per collision while 
b c is a dimensionless collisional average 

fe c = T7-^r [ dedld26{v 12 ■ g ) \v 12 ■ ef mim2 /(Pi,tn,t)/(p a ,g*,t), (10) 

2i0E c J mi + 7712 
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where it is understood that rrii — a i . The parameter b c is a dissipation parameter: 
7o& c is the ratio of the kinetic energy dissipated in an "average" collision to the 
instantaneous mean kinetic energy per particle. Eq. Q expresses the fact that due 
to the collisions the overall number of particles in the system increases in time. In 
some cases, this increase is so rapid that the gas exhibits a finite-time singularity, i. e. 
the number of grains become infinite in a finite amount of time (see next section). 
The existence of such a singularity, which can be thought of as a limitation of the 
model, signaling that some of the assumptions are too simple, or that some physical 
ingredient is missing, will depend on the parameters that characterize the evolution 
of the system, p, b c , a and the dimensionality D. 



3.1 Scaling behaviour 

We will analyze the evolution of the free evolving gas once the transient effects 
induced by the particular initial condition chosen have vanished. In particular, we 
will concentrate in the late time regime, addressing the issue whether the existence 
of a scaling solution is possible, where there exists a unique typical velocity and 
size of the particles characterizing the overall evolution of the system. 

To this end, a good way to characterize the possible kinetic scenarios, is to 
introduce an intrinsic time scale, C, which counts the number of collisions J7]. This 
scale is related to time through, dC = u)(t)dt. In terms of this variable, the evolution 
equations J7J and @ read 

dn , . 

^| = -2b c70 E c (II) 

with solutions 

n(t) = n e pC , E c (t) = E c0 < 3 - 2b " laC . (12) 

Assuming the existence of a scaling regime, the kinetics of the system is deter- 
mined by one single characteristic velocity v(t), and grain size f(t). The time depen- 
dence of the distribution function f(v,a,t) consequently occurs through these two 
quantities [up to a trivial time dependence through the density resulting from the 
normalization constraint 0], and the scaling assumption amounts to restricting / 
to the family of functions of the form: 

n ~ 

f(v,a,t) = f(c,a), (13) 

where we have introduced the rescaled size and velocity 

c » . 9= * . (14) 
v(t) f(t) 

A natural choice is to define v as the variance of the velocity distribution and f as 
the mean radius: v 2 — (v 2 ), f = (a). 

Within the scaling ansatz, we obtain from I1U1 that b c is time independent, and 
from (|SJ| that uj scales with time as nvf ' 1 (i.e. as the ratio of typical velocity over 
mean free path). If we use the fact that mass density is conserved (which implies 
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that r D n oc M oc t°), we know that r oc l/n 1//D and v 2 oc E c /M oc E c . Combining 
these relations, we get for the collision frequency 

- --D-1 ^ 1/0^1/2 f 1C N 

co oc vnr oc — oc n ' ii„' (15) 
r 

from which we can obtain a differential equation for the overall number of collisions 
as a function of time 




g-=^ exp[(g-& c7 o)c] (16) 



From this relation, we can obtain the explicit time dependence of the number of 
collisions, 

C W = -^_ln[l + (-g + 6c70 )^] (17) 

Finally, the explicit time dependence of the mean number of particles and kinetic 
energy read 



n(t) =n |l+(--^ + &c7o)^r x " ( 18 ) 



p/i-jy+bclo) 

2i>c7o/(-^+6c7o) 



E c {t) = E c0 \l + (-£ + b c7 o) ^otj (19) 

These expressions indicate when a finite-time singularity will develop. Due to mass 
conservation, the divergence in the number of particles implies that the mean size 
of the particles vanishes. We will see in the following sections that this feature may 
be associated with a "shattering" transition. In fact, whenever the inequality 

-|+Mo<0 (20) 
holds, the number of particles diverges at a time 

*/ = r P 1 , 1 (21) 

LOO — O c jo\ 

at which the kinetic energy vanishes. This inequality holds at large enough fragmen- 
tation probability, p > Db c yo- It therefore appears that the finite-time singularity 
cannot happen if fragmentation is just a marginal event (low p); alternatively, it 
will always be present in the elastic limit (a = 1 so that 70 = 0) . Hence, the kinetics 
of a gas of freely evolving fragmenting grains is determined by the inelasticity and 
the collisional average b c . In principle, this collisional average is itself a function of 
the inelasticity; however, as discussed in the Appendix, such a dependence only ap- 
pears when the velocity dependence of the distribution function f(v, a, t) deviates 
from a Maxwellian. 

In the absence of the above finite-time singularity, both n and E c show power 
laws with time, which, from (J7J in turn imposes that w k 1/t. The physical meaning 
of this remarkable feature (also observed in ballistic annihilation |5]) is that there 
exists a single relevant time scale in the problem. It is also noteworthy that from 
1151 . the constraint uj oc 1/t imposes a relation between the scaling exponents of n 
and E c . 



Fragmenting Granular Gases 



7 



The previous results presume the existence of a scaling regime, hypothesis that 
we have validated numerically, as will be described in Sect. 0] However, note that 
the existence of a finite-time singularity does not rule out the possibility of a scaling 
regime. In order to analyze quantitatively such a regime, it is useful to introduce a 
slower time scale, r 

dr = n(t)dt (22) 

In this "fake" time scale, the finite-time singularity is no longer present. This feature 
favours the numerical analysis of the asymptotic scaling behaviour of the mixture. 
The internal time scale is related to r through AC = Lodr/n. As a result, in terms 
of this new time scale, the number of collisions is expressed as 



C = 



bojo + (1 - s) P 



\n! 1 



bcjo + ( 1 - ) P 



LO T 



(23) 



Hence, using the exponential growth of both the number of particles and kinetic 
energy in terms of the number of collisions, we get 



n(r) = no |1 + 
E c (t) = E c0 { 1 + 



bcjo + 
b c fo 



D 



P/{bo~to + {l--h)p] 

.c,r> (24) 



LOqT 



-2D/ £>+i 



(25) 



When compared with l|18^ and l|19jl . one can see that, although qualitatively similar, 
the additional factors that characterize the exponents in this time scale ensure the 
absence of a finite-time singularity. 

The time evolution of the physical quantities of interest in the scaling regime, as 
depicted in 12 H and 12 "it . depend on the inelasticity coefficient, a, the fragmentation 
probability, p, the dimensionality, D, and the collisional average b c . From all of 
these, only b c is not known a priori. In the Appendix, we show how it is related to 
the particle distribution. In particular, it is a constant for a Maxwellian distribution 
of velocities with a single mass-independent temperature. In this limiting case, 
analytic expressions for the exponents can be obtained. However, as we will see, 
such a factorization is too simplistic in generic situations. The collisional average 
will therefore have to be determined numerically, and in all cases it is the only input 
required to predict the corresponding exponents governing the scaling behaviour of 
the freely evolving gas. 

Incidentally, the existence of a scaling regime can be used at one's advantage 
in DSMC. As has been explained previously, one of the numerical difficulties in 
studying the fragmentation of a freely evolving granular gas lies in the fact that 
the collision frequency decreases with time, making it harder to work with larger 
systems sizes (which is what inevitably happens in this case due to he fragmenta- 
tion) . On the contrary, working with the rescaled velocities c instead of the original 
v significantly speeds up the program, and allows to probe longer time scales i8;. 



3.2 Kinetic scenarios 

As an illustration of the various kinetic scenarios covered by the previous analysis, 
let us consider in more detail a few limiting cases related to the kinetics of freely 
evolving gases. 
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In the absence of fragmentation, p = 0, obviously the number of particles re- 
mains constant, the system corresponds to that of a monodisperse granular gas, 
and the energy decays as 

[1 + bc'-foujat] 2 

recovering Haff 's law [9J. Here, the collisional parameter b c is close to 1, with devia- 
tions originating from the non-Gaussian features of the velocity distribution |15) . In 
this particular case, the two time scales r and t are just proportional to each other, 
yielding the same exponent. The assumption of spatial homogeneity prevents us 
from addressing the late stage deviations from Haff 's law that have been described 
in freely evolving granular gases UJ. 

For an elastic gas, a = 1, the kinetic energy is conserved. If grains can still 
fragment with a probability p, the number of particles will always diverge at a 
finite time, 

n(t) = ^ =. (27) 

In this case it is easier to appreciate the advantage of introducing the slower time 
scale t. In terms of this variable, the number of particles will increase as 



i(t) = n 



1 - — ) /)~'u 7 



D/C-D-l) 



(28) 



In order to analyze the existence of an asymptotic regime (hypothesis upon which 
both 1271 and 12811 have been derived), it is more convenient to use the latter time 
scale, because it does not lead to a finite-time divergence. In this way, it is possible to 
make use of the standard techniques to address the appearance of power-law kinetic 
regimes, since the indefinite growth ensures the possibility to span several decades. 
Hence, once the algebraic regime is achieved, accurate values of the exponents can 
be derived as well. If such exponents are characterized, it is always possible to 
obtain their counterparts controlling the evolutions in the real time scale t. 

The technique of resorting to a slower time scale will always be useful for the 
present class of fragmentation problems, because the exponent appearing in 12411 . 
characterizing the increase in the grain number density, is always positive for any 
dimension larger than 1. For D — 1, the elastic limit is singular. In this case, the 
scale C is proportional to r and the number of particles increases exponentially 
in r rather than algebraically. This behavior is recovered, by performing the limit 
D —* 1 in 124H . and is due to the fact that in one dimension the collision frequency 
is proportional to the number of particles for our elastic mixture. 



4 Numerical simulations 
4.1 The method 

We have implemented the Direct Simulation Monte Carlo technique (DSMC), which 
offers a particle method to solve numerically a Boltzmann-type equation, such as 
This method was originally introduced to study rarefied gases |1U) . and has been 
extended to study granular gases 11 . It has been used in a number of different 
situations and it is easy to adapt to the present system, where the number of grains 
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increases with time. There exist already a number of applications to situations 
where the number of particles is not a conserved quantity, as it is the case in 
ballistic annihilation 2 . 

The assumption of spatial homogeneity simplifies further the implementation 
of DSMC, because it is enough to sample uniformly the relative orientation of the 
colliding pair, e, without the need to keep track of the positions of the particles. 
Sequentially, a pair of particles is selected at random. Also the vector e is chosen 
randomly on the unit sphere and the kinematic constraint is checked. In case the 
particles can collide, the cross section is computed and the collision takes place 
proportionally to its value. The decrease in velocity of the particles due to the 
inelasticity reduces the cross section, hence the simulations become more expensive 
as time evolves. After the collision is executed (according to 0), one of the two 
particles is chosen randomly, and it is broken with probability p. 

The collision frequency of each collision i-j can be computed on the basis of the 
cross-section, and hence the time is advanced accordingly at each collision 

^ + — (29) 

where n is the instantaneous number density. It is then possible to keep track of 
the kinetic evolution of the system, as well as gather statistics on the distribution 
of sizes and velocities. In the following section we will analyze both the kinetics and 
the - self-similar in time- polydisperse distributions that the free evolution of the 
gas gives rise to. 



4.2 Numerical results: Kinetics 

We have first analyzed the kinetics in the homogeneous regime to test the scaling 
hypothesis discussed in Sect. |3] To this end, we have used DSMC simulations to 
study the evolution of a two-dimensional granular gas of circular disks with constant 
restitution coefficient a. In all cases we start from an initial condition where all 
particles have the same size and their velocities are distributed according to a 
Maxwellian (monodisperse gas of elastic disks). As it has already been described 
in Sect. [5] we restrict ourselves to the simplest model where colliding particles 
fragment with a uniform probability independent of the collision properties, and 
where only one of the post-colliding disks fragments in two out-coming grains. 

In order to validate the existence of a scaling regime, we have computed the total 
number of particles, n(t), the collision frequency, uj(t), the mean velocity, (v 2 ) 1 ^ 2 , 
as well as the mean radius, (a) — f. We have also computed different moments of 
the size, (a n )(t), and velocity, (v 2 ) n ' 2 (t), distribution functions. Within the scaling 
picture, we expect (u™) cx (f) n and (\v\ p ) oc (v) p , as far as the time dependence 
is concerned. The possibility of a multi-scaling behaviour has been reviewed by 
Ben-Nairn and Krapivsky |12| . 

In Figure^, we display the time evolution of a few relevant quantities in the real 
time scale, t. For the parameters chosen, the system shows a finite-time singularity. 
As already indicated, such a phenomenon makes it difficult to assess the existence 
of a scaling regime, and complicates its analysis. On the contrary, when displayed 
in terms of the slower time scale r, as shown in Fig. 0d, one can see the existence 
of a scaling regime, which can be characterized. An alternative possibility to assess 
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the existence of a scaling regime without introducing r would be to check if (o p ) 
scales like {a p ) n/p . 

In Figure |3i we display different moments of the velocity distribution, as well 
as the mean collisional velocity (u co ii), defined as the mean velocity modulus of 
colliding partners, i.e restricting to those particles that are in precollisional config- 
urations (since a typical collision involves a particle that is "hotter" than the mean 
background, this quantity is larger than (v), as can be observed on the figure). Due 
to the inelasticity of the collisions, the typical velocity of the particles decays as 
a function of time. The figures show that such a decay is algebraic in r, and also 
that there is a single exponent characterizing different moments of the velocity. 
The same holds for the characteristic size of the particles as a function of time, 
as depicted in Figl^J. In this case however, the scaling does not hold for negative 
powers of the size. We will see in the next section how such a deviation from simple 
scaling relates to the peculiarities of the size distribution, but it does not imply a 
breaking of the scaling hypothesis; it only signals the fact that moments of negative 
order are controlled by the smallest species in the mixture. 

We can characterize the asymptotic regime by a set of exponents, that in the 
terms of r can be written generically as 

n(r) ~ t e 
E{t) ~ r" 27 ' 

<t(t) ~ r"' 3 '. (30) 

Once these exponents are known, their counterparts associated with the t scale 
follow. For example, using that dr = no(l + ar) e dt, for e' / 1 we have 



1 + ar 
1 + aro 



, n a(t-to) 
1 - e ) 



(31) 



from which we can derive the time evolution of any physical quantity. 

In Table we give the values of the exponents fitted numerically for the 2- 
dimensional gas at different values of the inelasticty parameter. We can determine 
the different exponents independently, but they are related to each other since all 
of them are eventually functions of the collisional average b c . It is in fact rather 
straightforward to check that, once one of the three exponents is known, the re- 
maining ones can be determined using the relationship between them as described 
in Sect. [3] 

In Table |5] we list the values of the collisional average b c computed from the 
numerical simulations. We also display the values of the exponent e' derived using 
the corresponding computed values of b c ; comparing such values with those mea- 
sured numerically (as listed in Table shows a good agreement, illustrating the 
consistency of the analysis. The scaling hypothesis also implies, from 112H . that 

f n \ 2b "~ m/p 

which provides an alternative way to compute b c in the scaling regime (where it 
has to be time independent). We display in Table |5] the collisional average com- 
puted following this route [under the symbol be ]■ It shows a good agreement with 
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Table 1. Exponents characterizing the time evolution of the number and energy 
densities as well as the mean particle radius (see (!■'->( )li for definitions) for different 
values of the inelasticity parameter a, when the probability of grain fragmentation 
at a collision is p = 1/2 



a 


e' 




i 


0.3 


0.95 


0.45 


0.55 


0.5 


1.05 


0.5 


0.475 


0.7 


1.3 


0.65 


0.35 


0.9 


1.9 


0.95 


0.125 


0.95 


2.3 


1.2 


0.05 



the computed collisional average, giving additional support to the existence of an 
asymptotic scaling regime. 

In the Appendix, we have given a generic expression for 6 C , showing that it 
depends sensitively on both the size and the velocity distribution of the mixture in 
the scaling regime. If the velocity dependence enters only through a Maxwellian, 
then the value of 6 C is independent of the mass distribution. This surprising result 
shows that this collisional average is a quantity that will be sensitive to deviations 
from Maxwell distributions. In this contribution we have only analyzed the behavior 
of the gas within the molecular chaos assumption, and hence, deviations of b c from 
unity signal a deviation from Maxwellian behavior (with the strong assumption that 
the temperature of a given species does not depend on its size). In a more generic 
study, deviations from Maxwellian behavior may also result from the breakdown of 
molecular chaos, so that the use of b c has to be supplemented by a second quantity 
that is sensitive to such a breakdown (e.g. the averaged impact parameter |13I14| 1. 

For a mono-disperse granular gas, deviations from Maxwellian behavior vanish 
as the elastic limit is approached and can be described quite accurately in terms of a 
Soninc expansion around the Maxwellian |15|. The values of b c displayed in Tabled 
on the contrary, show that the deviations from a Maxwellian velocity distribution 
become more pronounced as one approaches the elastic limit. This counter-intuitive 
trend can be traced back to the peculiar asymptotic size distribution of particles 
close to elasticity. Such a behavior is described in detail in the next section where 
it appears that the velocity distribution for nearly elastic gases deviates strongly 
from a global Maxwellian. At higher inelasticities, such deviations, even if always 
present, are not so dramatic. 

4.3 Numerical results: Distribution functions 

The existence of a scaling regime implies that the distribution function of a given 
quantity at two different times has to be invariant after an appropriate rescaling 
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Table 2. Collisional average b c , and predicted exponent e' for different values of the 
inelasticity parameter a, when the breaking probability is p = 1/2. be is obtained 
from l|32jl . plotting the kinetic energy as a function of the number of particles. 



Q 


b c 




4h 


0.3 


1.24 


1.248 


0.94 


0.5 


1.18 


1.173 


1.06 


0.7 


1.07 


1.042 


1.29 


0.9 


0.67 


0.645 


1.77 


0.95 


0.39 


0.2318 


1.94 



of this quantity [see (1131 and H14I 1. For instance, the distribution of a = cr/{a) 
is expected to be time independent. From the DSMC simulation runs, we have 
computed scaled distribution functions, which do indeed reach a steady state, in 
agreement with the scenario developed in Sect.|H] We will now discuss in detail the 
characteristic features of some of the relevant stationary distributions. 

Size distribution The fragmentation process generates a polydisperse mixture of 
grains that evolves continuously in time, with a decreasing mean size. Large particles 
are consequently progressively destroyed, but the shape of the size distribution is 
preserved, in such a way that the relative amount of particle with respect to the 
time dependent mean size, (cr)(£), becomes independent of time. 

In this regime, the scaled size distribution exhibits a marked peak for vanishingly 
small grains, indicating that smaller particles become predominant in the mixture. 
In Figure[!|]we show such distributions for a number of inelasticity values, where two 
regions can be identified. On the extreme of large grains, the decay is exponential 
while it is algebraic for small particles. The exponent characterizing this algebraic 
behavior increases with a. For example, for a — 0.3 the exponent is consistent with 
P(a) ~ a~° ' l , while for a — 0.7 we can fit the curve numerically with P(a) ~ cr -0 ' 6 . 

It is interesting to note that at a — 0.9, and also at a — 0.95, the above exponent 
becomes close to —1. Such a value is peculiar, because it implies that the probability 
distribution is not integrable at the origin, incompatible with the requirement that 
P{cr) [or P(<r)] has to be normalizable. In fact, a careful inspection of Fig. |3]shows 
that P{a) at high inelasticity develops a jump very close to the origin. Such a jump 
is consistent with the appearance of a singular contribution, due to a finite fraction 
of particles with vanishing size. This behaviour is reminiscent of a "shattering" 
transition £Q and is corroborated by the velocity distributions described in the 
next section. It indicates the development of an intrinsic size heterogeneity in the 
system. 

The algebraic decay of P(cr) implies that, even when the exponent is smaller 
(in absolute value) than —1, certain averages of the distribution function will not 
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exist. This is consistent with the behavior described in the previous section, and 
simply indicates that in general, averages of negative powers of the grain size will be 
controlled by the smallest species available in the mixture, rather than by the overall 
mixture. However, such a lack of scaling for certain moments of the distribution is 
not incompatible with the dynamic scaling behavior of the complete distribution 
function; rather, it shows some of the peculiarities of the asymptotic behavior of 
the scaling distribution (in this case at small sizes). 

The existence of a stationary size distribution suggests the possibility to map 
the fragmenting granular gas onto an effective polydisperse mixture with a time 
dependent mean. Although this is indeed possible in principle, we do not have 
any theoretical framework (except the complicated initial Boltzmann equation) to 
predict either the algebraic decay of P(a) or the details of the exponential tail. 

Velocity distribution We have also studied the velocity distribution to assess 
how relevant deviations from Maxwellian are, and also to analyze the meaning of 
a global temperature characterizing the whole system. 

In Figure 01 we display the distribution of rescaled velocities c at different in- 
elasticities. Deviations from Maxwellian behavior are observed at all values of a. 
However, the deviations become more prominent close to the elastic limit, in agree- 
ment with the development of a singular contribution coming from vanishingly 
small grains, explaining the small values for b c listed in Tabled In order to gain 
more insight, we have also computed separately the velocity distribution from those 
particles with a mass smaller than 5% the mean mass. In Figure we compare 
the overall velocity distribution with that of the smaller particles and the contri- 
bution coming from the rest of the grains for a = 0.95. The two subsets behave 
more closely to a Maxwellian than the mean velocity distribution. The low velocity 
part of the distribution is controlled by the contribution from most of the parti- 
cles, while the large velocity part is dominated by the small particles. These are 
characterized by a much larger temperature than the mean, while the rest of the 
species are much colder (large particle have suffered a low number of collisions (and 
hence of fragmentations), which requires that they were among the very "cold" ones 
initially) . 

Although less dramatic, the same trend is observed for all values of a. For 
example, in FigQjJj we display the equivalent distribution functions at a — 0.5. In 
this case, again the two subsets of grains evolve following separate approximate 
Maxwellians. According to this interpretation, as a increases, the separation in 
temperates between the two sets of particles increases. 

Such an inhomogeneous dependence of the velocity distribution (and in partic- 
ular its second moment) on grain species, explains the large deviations of b c from 
the Maxwellian prediction reported close to the elastic limit. These results also cast 
doubts on the use of a single temperature to characterize a polydisperse granular 
mixture. Nonetheless, it is worth noting that such a size heterogeneity does not 
affect the scaling hypothesis. Even if different temperatures could be defined for 
different subsets of particles, all of them have the same time evolution, explain- 
ing why it is possible to have a simple velocity scaling. Such a behavior has been 
reported for the free evolution of bidisperse granular mixtures |16|. However, the 
absence of a common temperature invalidates the use of perturbative calculations, 
e.g. through a Sonine expansion, to determine the deviations from Gaussian behav- 
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ior as the inelasticity is increased, as is customary in the case of monodisperse (and 
even bidisperse) granular gases. 

5 Discussion 

We have investigated the evolution of a granular gas in the absence of any driving, 
taking into account the possibility that the collisions may lead to a fragmentation 
of the impacting grains. 

Rather than performing a detailed analysis of the fragmentation process, we 
have analyzed the simplest possible model, where the fragmentation is a random 
event uncorrelated to the details of the collision, except for the requirement that 
only colliding particles may fragment. This simple fact already leads to a number of 
interesting results. We have shown the existence of a scaling homogeneous regime 
characterizing the kinetics of the system. Such a regime is controlled eventually by a 
single collision average, b c . The relevance of b c in the asymptotic regime shows that 
a theory based on rate equations, usually disregarding collision induced correlations, 
will be too simplistic to capture the dynamics in these systems. 

We have shown how the evolution in the scaling regime may lead to a finite-time 
singularity where the number of particles diverges. We have derived the conditions 
leading to this finite-time catastrophe, depending on the different parameters that 
characterize the system. In the absence of such a singularity, the various moments 
of the size and velocity distribution exhibit algebraic scaling laws with time, that 
can be regarded as a generalized Haff's law. 

The size distribution exhibits two different limiting kinds of behaviour: an ex- 
ponential tail at large sizes (compared to the time dependent mean size), and a 
power-law distribution for small sizes, which may develop into delta singularity 
close to the elastic limit. The details and character of such a transition, which has 
been reported in other fragmenting systems deserve a more careful analysis. 
The existence of a scaling regime is not incompatible with the presence of marked 
size heterogeneities in the mixture. In particular, the characteristic temperature of 
small and large grains can be significantly different, in such a way that the mean 
scaled velocity distribution cannot be regarded as a slightly distorted Maxwellian. 

In order to assess to which extent the results discussed are specific of the frag- 
mentation mechanism that has been chosen, we have also carried out numerical 
studies for a model where the fragmentation probability is correlated to the energy 
lost in a collision |17| . In that case there exists still a scaling regime characteriz- 
ing the kinetics. The scaled distribution functions are also similar, except for the 
fact that size distributions are closer to an exponential. Hence, the results obtained 
through the analysis of the simplest model provide a reference system, against which 
other models can be tested, and in this way it should be possible to disentangle the 
generic aspects involved in the dynamics of fragmenting granular materials. 
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associe" , and the hospitality of the "Laboratoire de Physique Theorique" in Orsay 
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Appendix: On the collisional average b c 

The kinetics of the free fragmenting granular gas depends on the value of the 
collisional average b c , which in turn is a complicated function of the joint size and 
velocity distribution. To proceed further, we assume that in the scaling regime, 
the velocity dependence is Maxwellian. From the distribution function f(v,a,t), 
we change variables and introduce the joint distribution of masses and velocity 
f(m,v,t) that we approximate by 

f(m,i, Vi,t) ~ ij)(rm,t)<j>(miVi,t) (33) 

where ip is the mass distribution, while (fiirmvf) is the Maxwellian, 

/ „ \ D/2 

4>(rn i v i ,t)=\—\ e « . (34) 

The time dependence of <j> is encoded in f3 — l/(feT), see @. The key assumption 
here is that ksT does not depend on the mass m, as would be the case in an 
equilibrium polydisperse system. We have 

b _ dld2de9(v!2 -e)\vi2 a^ 2 ~ 1 [L 12 ^(mi)(p(miVi)^(Tn2)(l>{rn2V2) _ 3g . 
J dl d2 Ae6(vi2 ■ e)o^ 2 ~ 1 ^l}{mi)4>{rniv\)^j){m2)4>{rn2V2) 

where dl stands for dui dmi. We can express the collisional average b c in dimen- 
sionless quantities, and also transform the velocities to express them in terms of 
the center-of-mass and relative coordinates. One arrives at 

2(3 /dcdwidm 2 c 3 a 1 p 2 -V 12 V>(m 1 )^(m 2 )(m 1 m 2 ) g / 2 e-^i2 c2 / 2 
c ~ (D + 1) / dcdmidmaccrg- V(m 1 )i/'(m 2 )(m 1 m 2 ) D / 2 e-^i2 c2 / 2 

P + 3 

Jdmidm2 (p^) 2 ^12 °"r2 _1 ^{rn\)ip(m2){mim,2) D /2 
j Am 1 Am,2 2 (J i2~ 1 4>("m 1 )-ip(m,2)(mim2) D / 2 

where use has been made of 

dse(c-s)\c-e]"^ r ^ '- e" (37) 

So, 6 e = 1 if the velocity distribution follows a Maxwellian. Hence, this collisional 
average is a good measure of deviations from Maxwellian behavior in polydisperse 
granular gases. 
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Fig. 1. Time evolution in linear- log scale of the total number of particles, mean size 
and characteristic velocities for a = 0.5, and a fragmentation probability p = 1/2 in 
D — 2. (a) Time evolution as a function of the real time scale t; (b) Time evolution 
in a log-log scale in terms of the slow time scale r. The thin lines correspond to the 
fitted power-laws for the different quantities. 



Fig. 2. Time evolution of different moments of the (a) velocity distribution func- 
tion, and (b) size distribution function for the same parameters as Fig. In (a) 
the time evolution of the mean velocity of colliding grains is also depicted. In both 
cases the slow time scale is used, and the corresponding fitted power-law decays are 
displayed as thin lines. 



Fig. 3. Log-log plot of the rescaled particle radius distribution, for a uniform ran- 
dom fragmentation probability p=l/2inD = 2. 



Fig. 4. Linear- log plot of the rescaled velocity distribution c = v/v. 



Fig. 5. Linear-log plot of particle velocity distributions, scaled by the mean ve- 
locity, for a uniform fragmentation probability. "Small" refers to the distribution 
of velocities for particles having a mass smaller than 5% the mean mass. "Large" 
corresponds to the distribution of the rest of the grains, (a) a — 0.95; (b) a = 0.5 
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